*Load up data file
use "`base'\Data Files and Code That Produced Them\P09 Inpatient Analysis File", replace
 
*These coefficients come from "P11 Estimate population for weights.sas" and are the slope of a regression of population on age in years. So convert age in months in this file
*into years then multiply by slope to get weights. Add one so it is effectively proportion and deals with 0 issue. Stata will do the rescaling.  
gen pop_all_np =  1 - 0.00926*(months_21/12)
gen pop_m      =  1 - 0.01296*(months_21/12)
gen pop_f_np   =  1 - 0.00537*(months_21/12)
 
local sufs = "all all_np f f_np m"  
local vars  = "visit  illness injury_or_alc  alcohol_any inj_by_self inj_by_oth inj_accident   "

*Create fitted lines on either side of RD
foreach suf of local sufs {
   foreach var of local vars {  
     *Fit the left side polynomial; 
      reg `var'_`suf'_r age_c age_c_sq  if age_c >= -2 & age_c < 0 [aweight=pop_`suf'], robust
      predict `var'_`suf'_left        if age_c >= -2 & age_c <= 0
     *Fit the right side polynomial;
      reg `var'_`suf'_r age_c age_c_sq  if age_c > 0 & age_c < 2 [aweight=pop_`suf'], robust
      predict `var'_`suf'_right       if age_c >= 0 & age_c < 2
   }
}
  
  
keep if age_months >= 21 - 2  & age_months < 21 + 2   

sum alcohol_any_all_np_r inj_by_self_all_np_r inj_by_oth_all_np_r inj_accident_all_np_r

*Total split into two groups
gen dif = injury_or_alc_all_np_r - (alcohol_any_all_np_r + inj_by_self_all_np_r + inj_by_oth_all_np_r + inj_accident_all_np_r)
sum  dif
drop dif
tab inj_accident_all_np_r

#delimit ;
graph twoway  
 	        (scatter inj_accident_all_np_r age_months, mcolor(blue)  msymbol(O)  msize(.6) yaxis(1) 
			   yscale(titlegap(2)) 
			   ylabel(0(35)70, axis(1) nogrid)) 
			(line inj_accident_all_np_left age_months, lwidth(thin) lcolor(blue) yaxis(1)) 
            (line inj_accident_all_np_right age_months, lwidth(thin) lcolor(blue) yaxis(1)) 
			  
	        (scatter alcohol_any_all_np_r age_months, mcolor(black)  msymbol(D)  msize(.4) yaxis(1))
            (line alcohol_any_all_np_left age_months, lwidth(thin) lcolor(black) yaxis(1)) 
            (line alcohol_any_all_np_right age_months, lwidth(thin) lcolor(black) yaxis(1)) 

	        (scatter inj_by_self_all_np_r age_months, mcolor(red)  msymbol(Th)  msize(.8) yaxis(1)) 
            (line inj_by_self_all_np_left age_months, lwidth(thin) lcolor(red) yaxis(1)) 
            (line inj_by_self_all_np_right age_months, lwidth(thin) lcolor(red) yaxis(1)) 
			
	        (scatter inj_by_oth_all_np_r age_months, mcolor(green)  msymbol(Sh)  msize(.8) yaxis(1)) 
            (line inj_by_oth_all_np_left age_months, lwidth(thin) lcolor(green) yaxis(1)) 
            (line inj_by_oth_all_np_right age_months, lwidth(thin) lcolor(green) yaxis(1)) 			
             , 
			 title("Figure 4: Inpatient Admissions by Cause",size(medlarge)) 
			 xtitle("Age at Time of Inpatient Admission",size(medlarge))  
			 ytitle("Rate per 10,000 People",axis(1) size(medlarge)) 
			 legend(off) 
			 text(62 20 "Accidental Injury",color(blue) size(medlarge)) 
			 text( 22 20 "Alcohol Intoxication", color(black) size(medlarge))
			 text( 7 22 "Deliberately Self Inflicted Injury", color(red) size(medlarge)) 
			 text( 15 22 "Deliberately Injured by Other Person", color(green) size(medlarge)) 
             graphregion(style(none) color(gs16))  
 
      ;
#delimit cr
grexportpdf "`base'\Code for Figures and Tables in Paper\Figure 4.pdf"


